The here() package

https://here.r-lib.org/

The here package creates paths relative to the top-level directory. The package displays the top-level of the current project on load or any time you call here():

here()

[1] "D:/Teaching/431/2023-431/431-slides-2023"

The here() package

Jenny Bryan’s Ode to the here package

Jenny Bryan on Project-oriendted workflow

Using R projects and the here package

How can you avoid setwd() at the top of every script?

1. Organize each logical project into a folder on your computer.
2. Make sure the top-level folder advertises itself as such. If you use RStudio and/or Git, those both leave characteristic files that get the job done.
3. Use the here() function to build the path when you read or write a file. Create paths relative to the top-level directory.
4. Whenever you work on this project, launch the R process from the project’s top-level directory.

Turning values like 77 and 99 into NA

Suppose we have the following small data set, where 77 = “Refused” and 99 = “No Response” or some other term that we want to think of as “missing”.

var1 <- c(20, 22, 35, 19, 77, 99)
var2 <- c(1, 3, 4, 77, 6, 99)
var3 <- c("Yes", "No", 77, 99, "No", "Yes")
dat <- tibble(var1, var2, var3) |> mutate(var3 = factor(var3))

miss_var_summary(dat)

# A tibble: 3 × 3
  variable n_miss pct_miss
  <chr>     <int>    <dbl>
1 var1          0        0
2 var2          0        0
3 var3          0        0

How can we convince R the 77s and 99s are missing values?

Use replace_with_na() from naniar

dat1 <- dat |>
  replace_with_na(
    replace = list(var1 = c(77, 99), var2 = c(77, 99),
                   var3 = c(77, 99)))
miss_var_summary(dat1)

# A tibble: 3 × 3
  variable n_miss pct_miss
  <chr>     <int>    <dbl>
1 var1          2     33.3
2 var2          2     33.3
3 var3          2     33.3

More on replace_with_na() here

Replacing 77 and 99 with NA across all variables

dat2 <- dat |>
  replace_with_na_all(
    condition = ~.x %in% c(77, 99))
miss_var_summary(dat2)

# A tibble: 3 × 3
  variable n_miss pct_miss
  <chr>     <int>    <dbl>
1 var1          2     33.3
2 var2          2     33.3
3 var3          2     33.3

Other ways to extend replace_with_na() are described here

What if we have the opposite issue?

The replace_na() function from the tidyr package (details here) replaces an NA value with a specified value.

In that sense, it is the compliment to the replace_with_na() function.

Demo: Replacing NA with a value

df <- tibble(x = c(1, 2, NA), y = c("a", NA, "b"))
df

# A tibble: 3 × 2
      x y    
  <dbl> <chr>
1     1 a    
2     2 <NA> 
3    NA b    

df1 <- df |> replace_na(list(x = 0, y = "unknown"))
df1

# A tibble: 3 × 2
      x y      
  <dbl> <chr>  
1     1 a      
2     2 unknown
3     0 b      

More on replace_na() here

Paired Samples: Comparing Proportions

Suppose we wish to investigate the association between retirement status and heart disease. One concern might be the age of the subjects: an older person is both more likely to be retired and more likely to have heart disease. In one study, therefore, 127 victims of cardiac arrest were matched on a number of characteristics that included age with 127 healthy control subjects: retirement status was then ascertained for each subject, with the results shown on the next slide.

Data from 254 subjects

We might think to use our usual twobytwo() function, but each matched pair in this study provides two responses, one for the Healthy Control and one for the Cardiac Arrest.

• Pearson \(\chi^2\) test doesn’t take this pairing into account.
• Our analysis must take this pairing into account.

Data from 127 matched pairs

The relevant data are stored for you in the retire.csv file.

The retire tibble

retire <- read_csv(here("c24", "data", "retire.csv")) |>
  mutate(healthy_control = 
           fct_relevel(factor(healthy_control), "Retired", "Not Retired"),
         cardiac_arrest = 
           fct_relevel(factor(cardiac_arrest), "Retired", "Not Retired"))

retire |> tabyl(healthy_control, cardiac_arrest) |>
  adorn_totals(where = c("row", "col")) |> adorn_title()

                 cardiac_arrest                  
 healthy_control        Retired Not Retired Total
         Retired             27          12    39
     Not Retired             20          68    88
           Total             47          80   127

• Goal: Estimate the relative odds of being retired for healthy individuals versus those who have experienced cardiac arrest.

What are we testing?

\(H_0\): There are equal numbers of pairs in which the healthy control retired and the matched individual who had a cardiac arrest did not retire, and in which the healthy control did not retire but the matched individual who had a cardiac arrest did retire.

• or \(H_0\): There is no association between health status and retirement.

Concordant vs. Discordant Pairs

The concordant pairs - or the pairs when two retired people or two non-retired people are matched - provide no information for testing a null hypothesis about differences in retirement status.

Thus, we focus only on the discordant pairs - or the pairs where a person who retires is paired with a person who has not retired.

The McNemar Odds Ratio

The McNemar odds ratio is just the ratio of the two discordant pairs, here: 12/20 = 0.60, which indicates the association’s strength – the odds of retiring for healthy controls are estimated to be 0.6 times as high as the odds of retiring for those with cardiac arrest.

McNemar’s test

\(H_0\): No association between retirement status and cardiac arrest

table(retire$healthy_control, retire$cardiac_arrest) |>
  mcnemar.test() 

    McNemar's Chi-squared test with continuity correction

data:  table(retire$healthy_control, retire$cardiac_arrest)
McNemar's chi-squared = 1.5312, df = 1, p-value = 0.2159

Without continuity correction, and tidied…

table(retire$healthy_control, retire$cardiac_arrest) |>
  mcnemar.test(correct = FALSE) |> tidy() |> gt()

McNemar’s test via epibasix package (includes odds ratio)

The mcNemar() function comes from the epibasix package.

table(retire$healthy_control, retire$cardiac_arrest) |>
  mcNemar(alpha = 0.1)

Matched Pairs Analysis: McNemar's Chi^2 Statistic and Odds Ratio
 
McNemar's Chi^2 Statistic (corrected for continuity) = 1.531 which has a p-value of: 0.216
 
McNemar's Odds Ratio (b/c): 0.6
90% Confidence Limits for the OR are: [0.28, 1.134]

McNemar’s test is typically valid when the number of discordant pairs exceeds 30.

McNemar test via exact2x2() (also includes odds ratio)

• The exact2x2 package can also do this, using a slightly different approximation for the confidence interval.

exact2x2(retire$healthy_control, retire$cardiac_arrest, paired = TRUE,
         conf.int = TRUE, conf.level = 0.90)

    Exact McNemar test (with central confidence intervals)

data:  retire$healthy_control and retire$cardiac_arrest
b = 12, c = 20, p-value = 0.2153
alternative hypothesis: true odds ratio is not equal to 1
90 percent confidence interval:
 0.3031389 1.1534984
sample estimates:
odds ratio 
       0.6 

Example 1 Setup (1/2)

In a double-blind trial, patients with active rheumatoid arthritis will be randomly assigned to receive one of two therapy types: a cheaper one, or a pricier one.

The proposed primary analysis will estimate the difference in the proportion of participants who had a DAS28 of 3.2 or less at week 48. The DAS28 is a composite index of the number of swollen and tender joints, the erythrocyte sedimentation rate, and a visual-analogue scale of patient-reported disease activity.

Example 1 Setup (2/2)

The study’s power analysis established a sample size target of 225 completed enrollments in each therapy group, based on a two-sided 10% significance level, and a desire for 90% power, assuming the proportion of participants with a DAS28 of 3.2 or less at week 48 was 0.27 under the less effective therapy.

What value was used in the power calculation for the proportion of participants with DAS28 of 3.2 or less at week 48 for the more effective therapy? Round your response to two decimal places.

Example 1: What code do we need?

• Are we comparing means or proportions?
• Paired or independent samples?
• What do we want to estimate? Necessary sample size? Power? Something else?
• What information is available to us?

Example 1 Code

power.prop.test(n = 225, p1 = 0.27, sig.level = 0.1, power = .9)

     Two-sample comparison of proportions power calculation 

              n = 225
             p1 = 0.27
             p2 = 0.3996866
      sig.level = 0.1
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

Rounding to two decimal places, p2 was 0.40

Example 2 Setup (1/4)

An investigator plans to replicate part of a study of the gut hormone fragment peptide \(YY_{3-36}\) (PYY) which reduces appetite and food intake when infused into subjects of normal weight. The original study is found here.

In common with the adipocyte hormone leptin, PYY reduces food intake by modulating appetite circuits in the hypothalamus. However, in obesity there is a marked resistance to the action of leptin, which greatly limits its therapeutic effectiveness.

Example 2 Setup (2/4)

The investigator wants to know whether obese subjects are also resistant to the anorectic effects of PYY. She intends to perform a randomized, placebo-controlled, double-blind crossover study on healthy obese subjects (including an equal number of male subjects and female subjects), with each subject studied on two occasions one week apart.

The subjects will be screened by a dietitian who will assess their eating behavior with (several established scales).

Example 2 Setup (3/4)

On two consecutive Thursdays, we will measure caloric intake during a buffet lunch offered two hours after the infusion of that week’s exposure. In one of the weeks, the subject will receive an infusion of PYY, and in the other week (with the order of the weeks determined at random) the subject will receive a placebo. The number of calories consumed at each lunch is measured and then converted to an appetite rating. Our primary outcome is the difference between the appetite rating after PYY and the appetite rating after placebo.

Example 2 Setup (4/4)

A clinically meaningful difference, the investigator tells you, would be one in which these comparisons would differ by 30 or more points on the appetite rating scale comparing the two infusions, which is 60% of the anticipated standard deviation of these results. The investigator then asks whether a study which gathers a total of 64 observations will yield at least 90% power to detect an effect of this size using a 5% two-tailed significance level, and to meet all other requirements described above.

Example 2: What code do we need?

• Are we comparing means or proportions?
• Paired or independent samples?
• What do we want to estimate? Necessary sample size? Power? Something else?
• Desired significance level \(\alpha\) = 0.05 (two-tailed)
• Sample Size is 64 observations (or 32 pairs)
• Delta is ?
• Standard Deviation is ?

Example 2 Code

power.t.test(n = 32, delta = 30, sd = 50, sig.level = 0.05, type = "paired")

     Paired t test power calculation 

              n = 32
          delta = 30
             sd = 50
      sig.level = 0.05
          power = 0.9077895
    alternative = two.sided

NOTE: n is number of *pairs*, sd is std.dev. of *differences* within pairs

What is the proper conclusion? Do we have at least 90% power?

ChatGPT

https://chat.openai.com/

ChatGPT response March 2023

ChatGPT 3.5 response November 2023

ASA Statement (2023-08-04)